clear
set mem 1g
set more off
prog drop _all
capture log close

capture log close 

pause off

	local graph1_options "ms(T) mc(black) mfc(black) c(l) lp(shortdash) lc(black)  lwidth(1.2)"
	local graph2_options "ms(O) mc(green) mfc(green) c(l) lp(dash_dot) lc(green)   lwidth(1.2)"
	local graph3_options "ms(T) mc(blue) mfc(none) c(l)   lp(longdash) lc(blue)   lwidth(1.2)"
	local graph4_options "ms(D) mc(red) mfc(red) c(l)     lp(solid) lc(red)   lwidth(1.2)"
	local graph5_options "ms(D) mc(red) mfc(none) c(l) lp(dash) lc(red)"
	local graph6_options "ms(O) mc(green) mfc(none) c(l) lp(dash) lc(green)"

display "$inout"
display "$inout_label"
set scheme s1color

****************************************************************

global inout "_CF_EducPlus_3"

global  dir_basedata    "out_files" 
global  dir_cfdata   	"out_files$inout"  
global  dir_work_out   	"out_figs$inout"  
global  dir_fig_out     "../outputs"


use $dir_cfdata/sim_model_part1, clear

rename * cf_*
rename cf_id  id 
rename cf_age age
*rename cf_cat_theta* cat_theta* 

merge 1:1 id age using $dir_basedata/sim_model_part1
keep if _merge == 3


** define variables start here:: 
gen cat_thetac = thetac_decile
gen cat_thetan = thetan_decile
tabstat thetac, by(cat_thetac) stat(mean N)
tabstat thetan, by(cat_thetan) stat(mean N)
label var cat_thetac "Cognitive ability deciles"
label var cat_thetan "Non-cognitive ability deciles"


** control for initial education 
sum age 
sort id age
by id: egen initeduc = min(educ)


foreach var of varlist thetac thetan  {

	sum `var' 
	gen `var'_raw = `var'
	egen `var'_std = std(`var') 

** standardize cognitive ability and noncognitive ability for easy interpretation

*	replace `var' = `var'_std

}

label var thetac "Cognitive ability"
label var thetan "Noncognitive ability"
label var thetac_std "Cognitive ability (standardized)"
label var thetan_std "Non-cognitive ability (standardized)"
label var initeduc "Initial years of schooling"
label var parents_clg "Parents are 4-year college graduates"

*gen thetaX = thetac * thetan 
*label var thetan "Cognitive ability $\times$ Noncognitive ability (standardized)"


foreach var of varlist thetac thetan {
	
	sum `var', detail
	gen `var'_lomed = (`var' <= `r(p50)' )
	gen `var'_himed = (`var' > `r(p50)' )

}

gen theta_hihi = thetac_himed * thetan_himed
gen theta_hilo = thetac_himed * thetan_lomed
gen theta_lohi = thetac_lomed * thetan_himed
gen theta_lolo = thetac_lomed * thetan_lomed



label var age "Age"
label var educ "years of schooling"	
label var dq  "smoking"	
label var addiction "Years of smoking"	
label var qstock "Addictive capital (q)"
label var cf_enroll "Enrollment in the counterfactual"
label var cf_dq     "Smoking status in the counterfactual"

	
	sort id age
	foreach var of varlist addiction cf_addiction dq cf_dq {
	foreach a in 20 25 30 {
		gen `var'_a`a'_ = `var' if age == `a'
		by id: egen `var'_a`a' = max(`var'_a`a'_)
		drop `var'_a`a'_
	}
	}
	
	
	capture drop *everdq* 
	
	sort id age
	capture label drop everdq
	label define everdq 0 "never smoker"  1 "ever smoker"  
	
	foreach a in 20 25 30 {
	by id: gen everdq_a`a' = (addiction_a`a' + dq_a`a' > 0)  
	label var everdq_a`a' "Ever Smoker at age `a' (Benchmark Case)"
	label values everdq_a`a' everdq
	}
	

	foreach a in 20 25 30 {
	by id: gen cf_everdq_a`a' = (cf_addiction_a`a' + cf_dq_a`a' > 0)  
	label var cf_everdq_a`a' "Ever Smoker at age `a' (Counterfactual)"
	label values cf_everdq_a`a' everdq

	tab everdq_a`a' cf_everdq_a`a' , row
	}
	

	capture label drop trgroup
	label define trgroup 0 "never smoking"  1 "switches to not smoking" 2 "always smoking" 
	capture drop trgroup*
	
	gen trgroup = . 		
	sum age
	foreach a of numlist `r(min)'/`r(max)' {

		** never smokers (control group)
		replace trgroup = 0 if dq == 0 & age == `a'
		** switchers 
		replace trgroup = 1 if dq == 1 & cf_dq  == 0 & age == `a'		
		** always smokers
		replace trgroup = 2 if dq == 1 & cf_dq == 1	& age == `a'		
		
	}
	label values trgroup trgroup
	
	foreach var of varlist trgroup {
	foreach a of numlist 20/30 {
		qui gen `var'_a`a'_ = `var' if age == `a'
		qui by id: egen `var'_a`a' = max(`var'_a`a'_)
		qui drop `var'_a`a'_
	}
	}
	





***********************************************************************************

capture program drop plot_chg_BY_age 
program plot_chg_BY_age 
	
	args var 
	
   capture drop chg_`var' 
   capture drop pchg_`var'  
      

   local ly : variable label `var'

   gen chg_`var' = cf_`var' - `var' 

   gen pchg_`var' = chg_`var'/`var'
   

   label var chg_`var'  "Treatment effect on `ly'"

   label var pchg_`var' "Treatment effect on `ly' (%)"
	
	
	local graph1_options "lp(shortdash) lc(black) c(l) "
	local graph2_options "ms(*) mc(blue) mfc(blue) c(l) lp(dash_dot) lc(blue)"

*	local graph1_options "ms(T) mc(blue) mfc(blue) c(l) lp(solid) lc(blue)"
*	local graph2_options "ms(D) mc(red) mfc(red) c(l) lp(solid) lc(red)"
	local graph3_options "ms(O) mc(green) mfc(green) c(l) lp(solid) lc(green)"
	local graph4_options "ms(T) mc(blue) mfc(none) c(l) lp(dash) lc(blue)"
	local graph5_options "ms(D) mc(red) mfc(none) c(l) lp(dash) lc(red)"
	local graph6_options "ms(O) mc(green) mfc(none) c(l) lp(dash) lc(green)"

		
	** treatment by age 
	preserve 
	
	collapse *chg_`var' , by (age)

    label var chg_`var'  "Treatment effect on `ly'"
    label var pchg_`var'  "Treatment effect on `ly' (%)"
    	
   	graph twoway scatter chg_`var' age if age<=30, `graph2_options'  scheme(s1color) ylabel(, grid)
	graph export $dir_work_out/treatment_chg_`var'_BY_age30.pdf, as(pdf) replace 	
	graph export $dir_work_out/treatment_chg_`var'_BY_age30.eps, as(eps) replace 	
	

	restore 	
	
end 


capture program drop plot_chg_BY_theta 
program plot_chg_BY_theta  
	
	args var 

   capture drop chg_`var' 
   capture drop pchg_`var'  
   capture drop revchg_`var'
	
   local ly : variable label `var'
  
   gen chg_`var' = cf_`var' - `var' 
	
*   gen pchg_`var' = chg_`var'/`var'
   
   gen revchg_`var' = `var' - cf_`var' 
	
	local graph1_options "lp(shortdash) lc(black) c(l) "
	local graph2_options "ms(*) mc(blue) mfc(blue) c(l) lp(dash_dot) lc(blue)"

*	local graph1_options "ms(T) mc(blue) mfc(blue) c(l) lp(solid) lc(blue)"
*	local graph2_options "ms(D) mc(red) mfc(red) c(l) lp(solid) lc(red)"
	local graph3_options "ms(O) mc(green) mfc(green) c(l) lp(solid) lc(green)"
	local graph4_options "ms(T) mc(blue) mfc(none) c(l) lp(dash) lc(blue)"
	local graph5_options "ms(D) mc(red) mfc(none) c(l) lp(dash) lc(red)"
	local graph6_options "ms(O) mc(green) mfc(none) c(l) lp(dash) lc(green)"

	foreach figage of numlist 30 {
	
	preserve 
	
	keep if age == `figage'
	
	** This is ATE, not the ratio between median's 
	sum chg_`var'  
	local ate = `r(mean)'
	
	sum revchg_`var' 
	local revate = `r(mean)'
	
*	sum pchg_`var'  
*	local pate = `r(mean)'
	
	collapse *chg_`var', by (cat_thetac cat_thetan)	
	
	sum *chg_`var'
	
	label var chg_`var'  "Treatment effect on `ly'"
*	label var pchg_`var'  "Treatment effect on `ly' (%)"
	label var revchg_`var'  "Treatment effect on `ly'"	

*	_pctile chg_`var', nquantiles(5) 
*	twoway contour chg_`var' cat_thetac cat_thetan, ccut(`ate' `r(r1)' `r(r2)' `r(r3)' `r(r4)' `r(r5)') crule(linear) ecolor(green*0.1) scolor(green*1) zlabel(, format(%9.2f))			


	local maxn = 5 
	local maxnp = `maxn'+1
	
	
	sum chg_`var'
	local rmax = `r(max)'
	local rmin = `r(min)'	
	
	if (`ate' <= 0){
*		 local rmax = 0
	} 
	else {
*		local rmin = 0
	}
	
	_pctile chg_`var', nquantiles(`maxnp') 
	return list 

*	local rmax = r(r`maxn')
*	local rmin = r(r1)		

	local dis = 1000 
	local n_dis = 1
	
	foreach n of numlist 0/`maxn' {
		
		local r    = `rmin' + `n' * (`rmax' - `rmin')/`maxn'
		
	*	local r = r(r`n')
		
		local r`n' = `r'
		if ( abs( `r' - `ate') < `dis' ){
			local dis = abs( `r' - `ate')
			local n_dis = `n'
*			disp `dis'
*			disp `n_dis'
		}
		disp `r'
	}
	local r`n_dis' = `ate'
	
	if (`ate' <= 0){
		twoway contour chg_`var' cat_thetac cat_thetan, ccut(`r0' `r1' `r2' `r3' `r4' `r5' ) crule(linear) ecolor(green*0.1) scolor(green*1) zlabel(, format(%9.2f))	interp(none)	scheme(s1color) 	
	} 
	else { 
		local r0 = 0.0 
		twoway contour chg_`var' cat_thetac cat_thetan, ccut(`r0' `r1' `r2' `r3' `r4' `r5' ) crule(linear) scolor(purple*0.1) ecolor(purple*1) zlabel(, format(%9.2f)) 	interp(none)	scheme(s1color) 	
	}
	
	graph export $dir_work_out/treatment2_chg_`var'_BY_cat_theta_age`figage'.pdf, as(pdf) replace 
	graph export $dir_work_out/treatment2_chg_`var'_BY_cat_theta_age`figage'.eps, as(eps) replace 	
	
	restore 
	
	}
	
end 


***********************************************************************************



global yinf   -1
global ystep  0.2
global ysup   1

plot_chg_BY_age dq 
graph export $dir_fig_out/Figure_05a.pdf, as(pdf) replace 
graph export $dir_fig_out/Figure_05a.eps, as(eps) replace 


plot_chg_BY_theta dq 
graph export $dir_fig_out/Figure_05b.pdf, as(pdf) replace 
graph export $dir_fig_out/Figure_05b.eps, as(eps) replace


capture log close

